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Abstract 

We present algorithms and heuristics to compute the characteristic 
polynomial of a matrix given its minimal polynomial. The matrix is rep- 
resented as a black-box, i.e., by a function to compute its matrix-vector 
product. The methods apply to matrices either over the integers or over 
a large enough finite field. Experiments show that these methods perform 
efficiently in practice. Combined in an adaptive strategy, these algorithms 
reach significant speedups in practice for some integer matrices arising in 
an application from graph theory. 

Keywords: Characteristic polynomial ; black-box matrix ; finite field. 

1 Introduction 

Computing the characteristic polynomial of an integer matrix is a classical math- 
ematical problem. It is closely related to the computation of the Frobenius nor- 
mal form which can be used to test two matrices for similarity, or computing 
invariant subspaces under the action of the matrix. Although the Frobenius 
normal form contains more information on the matrix than the characteristic 
polynomial, most algorithms to compute it are based on computations of char- 
acteristic polynomials (see for example [25, §9.7]). 

Several matrix representations are used in computational linear algebra. In 
the dense representation, amxn matrix is considered as the array of all the mxn 
coefficients. The sparse representation only considers non-zero coefficients using 
different possible data structures. In the black-box representation, the matrix 
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is viewed as a linear operator, and no other operation than the application to a 
vector is allowed. Though constraining, this limitation preserves the structure 
or sparsity of the matrix and is therefore especially well suited for very large 
sparse or structured matrices. 

Computation of the characteristic polynomial of dense matrices has already 
been well studied both in theory and practice: over a finite field, [23, 24] set 
the best complexity (using respectively a deterministic and a probabilistic algo- 
rithm), and [8, 24] propose efficient implementations. Over the integers, the best 
complexity is achieved in [22], but currently the most efficient implementations 
rely on the Chinese remainder algorithm [8]. 

In the latter article, a competitive approach is introduced that limits the 
use of the Chinese remainder algorithm to the computation of the minimal 
polynomial. The characteristic polynomial is then recovered by determining the 
multiplicities of its irreducible factors. This task is done using the deterministic 
algorithm for the characteristic polynomial over a randomly chosen prime field. 

In the black-box model, the minimal polynomial is used as a building block 
for many algorithms over a finite field. Adapted from the iterative numerical 
methods (Lanczos, Krylov), the Wiedemann minimal polynomial algorithm [29, 
21] has excellent asymptotic complexity and is used in efficient black-box linear 
algebra software such as LinBox 1 . 

However less is known concerning the characteristic polynomial of black- 
box matrices. It is an open problem [20, Open Problem 3] to compute the 
characteristic polynomial as efficiently as the minimal polynomial, using the 
Wiedemann method. The latter uses 0(n) products of a square n x n matrix 
by a vector and O(n 2 (\ogn)°^) additional arithmetic operations with 0(n) 
extra memory storage. Eberly gives an algorithm using 0(n) matrix vector 
products, and (D(<f)n 2 ) additional operations, where (j> is the number of invariant 
factors of the matrix [12]. In the worst case, <f> — 0(n) and the algorithm does 
not improve on the complexity of dense algorithms. Villard proposed in [28] a 
black box algorithm to compute the Frobenius normal form and therefore the 
characteristic polynomial in 0{yfn log(n)) computations of minimal polynomials 
and C(n 2 5 (logn) 2 loglogn) additional field operations. 

Instead, we propose here several algorithms and heuristics focusing on effi- 
ciency in practice. The general strategy is to compute the minimal polynomial 
using Wiedemann's algorithm and decompose it into irreducible factors. There 
only remains to determine to which multiplicity each of these factors appear in 
the characteristic polynomial. In section 2 we propose several methods to deter- 
mine these multiplicities. Adaptive combination of them is discussed in section 
3. Under a conjectured hypothesis the latter is shown to require O {n^/n) ma- 
trix vector products which improves by a logarithmic factor on the complexity 
of Villard's algorithm. 

Lastly, an algorithm for the computation over the ring of integers is derived 
in section 4. It is based on the multifactor Hcnscl lifting of a gcd-free basis, 
following Storjohann [26] . The benefit of this approach is verified by experiments 
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presented in section 5. Several sparse matrices are considered, including a set 
of adjacency matrices of strongly regular graphs, coming from an application in 
graph theory. 



2 Three methods for computing multiplicities 

In this section we consider a matrix A over a finite field K = GF(q). Let 
-Fmin = n»=i ^T* b e * ne decomposition of the minimal polynomial of A in 
irreducible monic factors. The characteristic polynomial is then 

^ ar =n p r 

i=l 

for some mj > ej. We also denote by the degrees of the factors: di = deg(Pj). 

To recover the multiplicities mj, we will present three techniques, based on 
black-box computations with the matrix A: the nullity method (§2.1) uses the 
rank of Pi (A) to reveal information on the multiplicity mi, the combinatorial 
search (§2.2) is a branch and bound technique to solve the total degree equation 
whose integral unknowns are the multiplicities and the index calculus technique 
(§2.3) uses a linear system solving based on the discrete logarithm of equation 
1 evaluated in random values. 



2.1 The nullity method 

Deflniton 2.1. The nullity v{A) of a matrix A is the dimension of its nulls-pace. 



We also recall the following definitions: 

The companion matrix of the monic polynomial P = X d + J2i=o a iX l is 



the matrix C r 



o 
1 o 



-a 



Its minimal polynomial and its characteristic 



1 -a d -i 

polynomial are equal to P. 

The block Jordan matrix of an irreducible polynomial P of degree d to a 

"C P B 

power k is the kd x kd matrix J P k of the form J P k = 

Cp B 

L c P j 

where the d x d matrix B is filled with zeros except for Bd,\ = 1. Its minimal 
polynomial and its characteristic polynomial are equal to P k . This definition 
extends the usual notion of Jordan blocks for d = 1. 

The Frobenius normal form of a Matrix A is the unique block diagonal 
matrix F = Diag(C/ , Cj x , . . . ) such that A = U~ 1 FU for a nonsingular matrix 
U. The polynomials fi are the invariant factors of A and satisfy /o = P^in an d 
fi + i divides fi for all i > 0. 

The primary form of a Matrix A (also called the second Frobenius form 
in [15]) is a further decomposition of the Frobenius normal form where each 
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companion block C f i is replaced by a block diagonal matrix Diag( J^i , J^k 2 ,•■•)• 
The gj are the irreducible factors of ft, with the respective multiplicities kj. The 
primary form is unique up to the order of the blocks. 

Example 2.2. Consider the matrix in Frobenius normal form 

A = -Dia(/(C X 5_ 6X 4 + i4x3-i6x 2 +9x-2) C X 2_ 2 x+i) 
over GF(5). The corresponding primary form is the matrix 

B = Diag(J( X 2_2x-ip,Jx-2,J X 2-2x-i)- 



A = 





10 

10 

1 

1 



2 

-9 
16 
-14 
6 



r 
10 
1 
1 



The method of the nullity is based on the following lemma: 

Lemma 2.3. Let A be a square matrix and let P be an irreducible polynomial 
of degree d, of multiplicity e in the minimal polynomial of A, and of multiplicity 
m in the characteristic polynomial of A. Then v{P e (A)) = md. 

Proof. Let F be the primary form of A over K: F = U~ 1 AU for a non sin- 
gular matrix U. F is block diagonal of the form Diag(J p =j). Then P e (A) = 

U- 1 P e (F)U = C/- 1 Diag(P e (J plt ))C/. On one hand P e annihilates the blocks 

3 

Jpk where P = Pj and k < e. On the other hand, the rank of P e {J P k) is full for 

3 3 

P 7^ Pj, since P and Pj are relatively prime. Thus the nullity of P e {A) exactly 

corresponds to the total dimension of the blocks J P k where P 

j 

md. 



Pj, which is 



□ 

From this lemma the following algorithm, computing the multiplicity rrn of 
an irreducible factor Pi is straight-forward: 



Algorithm 1: Nullity 
Data: A: an n x n matrix over a Field K, 
Data: P: an irreducible factor of P^ ln , 
Data: e: the multiplicity of P in P^ in , 
Result: m: the multiplicity of P in P c ^ ar - 
begin 

r = rank(P e (A)) 

return m = (n — r) /degree(P) 
end 



Proposition 2.4. Algorithm 1 computes the multiplicity of P in the character- 
istic polynomial of an nx n matrix A using O(edn0,) field operations, where ft 
is the cost of a matrix-vector product with A, d is the degree of P and e is its 
multiplicity in the minimal polynomial. 
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Proof. Using Horner's rule, the matrix P(A) can be written as aoI n + A(a\I n + 
Ma-il n + ...)). Hence, applying a vector to this blackbox only requires d ap- 
plications of a vector to the blackbox A, ie. (D(dSl) field operations. Thus 
applying a vector to P e (A) costs 0(edSY). Lastly, the rank of this matrix can be 
computed in 0(ednfl) field operations, using Wiedemann's algorithm combined 
with preconditioners [11]. □ 

This algorithm is therefore suitable for irreducible factors P where the prod- 
uct ed is small. 

Now if e is large, the computation of rank(P e (^4)) may be too expensive. 
Still, some partial knowledge on the multiplicity can be recovered from the rank 
of the first powers of P(A). This can help to shorten the computation of other 
algorithms, as will be shown in section 3. We now describe how these partial 
multiplicities can be recovered. 

The multiplicity m of P is formed by the contribution of several blocks of 
the type J P j for j £ [1 . . . e] in the primary form of A. Whereas the blocks with 
small j can be numerous, there must be few blocks with large j, due to the 
limitation of the total dimension (since e is large). 

We denote by nij the number of occurrences of J p j in the primary form of 

i 

A. From the determination of the riij, we can directly deduce the multiplicity 
rrii by the relation 

e 

i=i 

We now show how to compute the riij for small j, using algorithm 1. 

Lemma 2.5. Let P be an irreducible polynomial of degree d over a finite field 
K and k and e > 1 be two integers. Then v{P k (Jp<:)) = min(k,e)d 

Proof. Let A = J P ^ and B = P k (A). If k > e, then P k is a multiple of the 
minimal polynomial of A. Thus B is the zero matrix, and its nullity equals its 
dimension: ed. 

Now suppose k < e. Let K be an extension of K such that P splits into 
d degree one factors Pi over K. Since any finite field is a perfect field, these 
factors are distinct. 

The minimal polynomial of A over K is still P e . Consequently, the Frobcnius 
normal form of A over K is Cpc and its primary form is F = Diag(Jpe). More 
precisely, there exists U <G M n (K) such that A = U~ 1 FU. We have therefore 
B = U- 1 P k {F)U = C/- 1 Diag(P fc (Jpe))C/. 

First consider the case k = 1: the minimal polynomial of each Pi(Jp^) is 
X e and so is the minimal polynomial of each P(Jp^) (since the Pi are relatively 
prime). Hence the primary form of P(Jp») is Jx=- Therefore there exist V € 
M n (K) such that 

B = U^V^DiagiJxe, J X «)VU. 

" v ' 

d times 
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Lastly the rank of Jx« being e — 1, we deduce that rank(B) = die — 1). The 
nullity of B is therefore v{B) = d. 
For the general case, we have 

B = U~ V- 1 Diag((J X e) fe , . . . , (J x ,) k )VU. 

V v ' 

d times 

Now Jx" is e x e and nilpotent with ones on the super-diagonal so that its fc-th 
power has rank max(0,e — k). Thus, rank(£>) = max(0, e — k)d and v(B) = 
min(e, k)d. □ 

We now apply this result to the irreducible factors of the minimal polynomial 
and denote the nullity of P? (A) by v itj = v{P((A)). 

First, the nullity of Pi(A), can be decomposed into the sum of the nullities 
of each Pi{J P k) for every k < ec 

= ^n itk di (3) 

k=l 

Now applying P? to A, every P^(J P k) for k < j will be a zero matrix and 
therefore contribute with kdi to the nullity. Otherwise, if k > j, the contribution 
to the nullity remains jd*. Therefore we have: 

3 e 4 

v%,i = y^^rij^kdj + rii^jdi (4) 

fe=l k=j + l 

From these two equations, we deduce the mj: first we have 

1 1 j ~ 1 
7^,j-i = 7 ^ n l k kdi + n t jdi + V] n^fed,. 

J J fe=i fe=j+l 

Now, since: YTk=j+i n i,kdi = fij+i — Vij, the number of occurrences directly 
is: 

1 ( 1 \ 1 , 

/J fc=1 

Therefore we obtain corollary 2.6 giving the expression of the riij: 
Corollary 2.6. 

ni,i = (2v iA - v ii2 )/di 
1 f 1 



1 i_1 

? — ^ Vj e [1 . . . ej] 



fe=i 

e 4 -l 



^ e t fc=i 



- ^2 ni > kk 
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The last formula for n^ ei is given for the sake of completeness: in practice, 
one will never compute every riij, since one would rather directly compute the 
nullity of Pf L l (A) instead, to deduce the multiplicity m t from algorithm 1 . 

2.2 The combinatorial search 

In the following, we want to determine the values of the unknown riij. They 
must satisfy the total degree equation: 



We can also discriminate potential candidates using the trace: the degree 
n — 1 coefficient of the characteristic polynomial is the negative of the trace of 
the matrix. Denote by ti the degree n — 1 coefficient of an irreducible factor 

P, = X di +t i X di - 1 + . . .. Then the degree n-1 coefficient of n^r* isE^ m *- 
We thus have the trace test: 



In a pure black-box model, the trace can be computed using n matrix-vector 
products. For many sparse or structured matrix representations, a faster method 
is available as well. 

Then it suffices to use e.g. a Branch-and-Cut algorithm to compute all the 
integer fc-tuples satisfying both equations (5) and (6). Of course, if some of the 
unknowns already computed (e.g. by the nullity method) the set of 

candidates is accordingly reduced. 

The remaining candidates will then be discriminated by evaluations of the 
characteristic polynomial at random values, i.e. computations of determinants 
of XI — A matrices. Indeed, we have efficient methods of computing the deter- 
minant of a black-box matrix (see e.g. [27, §3.1 Determinant Preserving Pre- 
conditioners] and references therein). Algorithm 2 sums up this combinatorial 
search strategy. 

2.3 Index calculus method 

Evaluating equation (1) at a point A leads to an equation over the finite field, 
where the multiplicities rm are the unknowns. Inspired by index calculus tech- 
niques [4] , the idea here is to consider the discrete logarithm of such an equation 
(with an arbitrary choice of generator), to produce a linear equation in the mj. 
Taking several of these equations for different Aj forms a linear system of equa- 
tions, with dimension k, the number of unknown multiplicities. 

In this discussion the base field is GF(q) and q is sufficiently large with 
respect to n as discussed below. The characteristic polynomial evaluated at a 




(5) 




(6) 
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Algorithm 2: Combinatorial-search 



Data: A, an n x n matrix 

Data: D = (di)i, the degrees of the irreducible factors Pi of P c h ar 
Data: M, a set of precomputcd riij 
Result: N = (riij) 
begin 

/* using degree and trace constraints */ 
sol = Branch-and-Cut(A, D, M) 
while #sol > 1 do 

Pick A € K at random 
5 = det(A7 - A) 
|_ Discard any TV € sol s.t. n( P /)" M ( A ) 7^ 6 
return N = sol[l] 
end 



given value A presents this equation in the unknown exponents mj : 

k 

f[pp(X)=det(XI-A). (7) 

Now, if A is not a root of the characteristic polynomial, taking the discrete 
logarithm of these terms for a generator g of the field leads to this equation 
modulo q — 1 : 

k 

^m J log g (P J (A)) = log g (dct(A/-A)) mod q - 1, (8) 

which is linear in the unknowns rrij. We can therefore build a I x k linear 
system by randomly choosing I values A^. This system is consistent since the 
multiplicities mi are a solution vector of this system. If the solution is unique, 
then it is the vector of multiplicities over Z. 

The computation of this vector can either be done by a dense Gaussian 
elimination over the ring Z 9 _i or over a finite field Z p where p is a large prime 
factor of q — 1 (larger than n). In this last case, the result will be correct as 
long as the system remains nonsingular modulo p. 

Algorithm 3 describes this techniques in more details. 

Let k be the number of unknown multiplicities. The first step is to find k 
values Aj forming a non singular system. Therefore, we propose, in algorithm 3, 
to evaluate the system at more than k points. The complexity of forming a row of 
G costs only 3 $^ j=1 di arithmetic operations using Horner's method. Therefore 
trying as many as n different values for A is a negligible cost. Furthermore, one 
could use fast multi-point evaluation to get blocks of rows simultaneously (up 
to n rows at a cost essentially linear in n). 
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Algorithm 3: Index-calculus 



Data: A, an n x n matrix over a finite field K = GF(q) 

Data: Pj, the irreducible factors of P^ in , 

Data: S, the set of indices of the unknown multiplicities, 

Data: Q = Ylj^s ' ^ e P ar ^ a l product of the irreducible factors with 

known multiplicity nij. 
Result: P c ^ ar or "fail" 
begin 

k = #S;l = 0;H=[]; 
Choose a generator g of K 
Let p be a prime factor of q 
while rank(H) < k do 

I = I + 1; if I > n then return "fail" 

Choose randomly A; e K 

repeat 

a lJ = Pj(^l) f° r au 3 

n = QW) 

until 7^0 and aj ^ 0, Vj e 5 

/* H is extended to the size Z x fc */ 

Stack the row [log g aij mod (q — 1)] mod p to H 

JC = { indices of the first k independent rows of H} 
Set B = a k x k nonsingular matrix of these rows 
Compute b = [log g (det(AjJ - A)) - log^)]^ 
Solve Bx = b 
return P = U je s P j i Q 
end 
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The rank of H is computed all along the process, each new row being in- 
crementally added to the triangular decomposition of the current matrix. This 
Gaussian elimination (performed by the LQUP algorithm [19] for example) also 
provides the indices of the first k linearly independent rows, and therefore the 
indices of the convenient A^. Lastly the vector b is formed, using only k deter- 
minant computations. 

In practice, it appeared, as in index calculus [2, 18], that the number of rows 
required to get a full rank matrix B is always quite close to k. However, we do 
not have a proof of this property, and we therefore state it as the conjecture 2.7. 

Conjecture 2.7. Let A be anxn matrix over Z and Pi ... Pj, be the irreducible 
factors of its minimal polynomial. Let p > n be a prime chosen randomly in 
finite set. Let q = 1 + \p of the form r k where r is a prime number. Let g be a 
generator of GF(q) and (Ai, . . . , A„) uniformly chosen at random in GF(q). Let 
H = [hij] where hij — (log fl (Pj(Aj)) mod (q — 1)) mod p. Then rank(H) = k 
with high probability. 

Informally, in our system the evaluations at the Aj are independent and can 
be considered as seeds for the pseudo-random generator of taking the discrete 
logarithm of the polynomial evaluation. Therefore the entries of the system 
modulo q—l are at least close to random entries as soon as the polynomials are 
distinct. Would they be true random values, the singularity /nonsingularity of 
the matrices would follow the analysis of e.g. [3, Corollary 2.4]: if L is a square 
matrix of uniformly random entries modulo q—l and p is a prime diving q—l, 
then the probability that L mod p is singular is of order ^ . 

Theorem 2.8. Assuming conj. 2.7, algorithm 3 is correct and its asymptotic 
complexity is 0{knVL) where VL > n is the cost of a multiplication of A by a 
vector. 

Proof. Let I be the number of rows required to get an invertible system, I > k. 
Each determinant computation requires 0(n) application of A to a vector [27, 
§3.1]. Building each row of the matrix requires a Horner like evaluation of a 
polynomial with total degree less than n, it therefore costs 0(n) operations. 
Triangularization of G requires C(/fc w_1 ) operations. Solving the system Bx = 
b, knowing the triangular decomposition of B requires 0(k 2 ) operations. The 
discrete logarithms can be tabulated [6] with 0{q) memory (or to avoid this 
extra memory, one can compute the whole sequence of powers of a generator of 
K , sort the matrix and vector entries and find the correspondences with some 
0~(lk + q) extra field operations). The overall complexity is thus O(kn0, + In + 
Ik^- 1 ) which is O(knn) when I = C(n^fc 2 -"). □ 

The algorithm stops arbitrarily when I = n + 1. We see here that a larger 
I is acceptable for the complexity result, but in our experiments a very small 
I (e.g. I « k) always suffices. In the following sections, this algorithm will be 
used with k < y/n, thus giving an expected 0{n 15 £l) complexity. 
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3 Adaptive black-box algorithm over a finite field 



We show in the present section how to combine the ideas of the previous section 
together with already existing techniques to form an adaptive algorithm com- 
puting the characteristic polynomial of a black-box matrix over a finite field. 
The algorithm is adaptive in the sense of [5], meaning that it chooses the best 
variant depending on discovered properties of its input. 

We first combine the nullity method with the combinatorial search. We 
then show an algorithm improving on the asymptotic complexity. Finally we 
give some improvements which are efficient in practice on typical matrices. 

3.1 Nullity method and combinatorial search 

These two algorithms are complementary: the nullity is efficient for the determi- 
nation of the multiplicites of factors of small degree, whereas the combinatorial 
search is adapted to the large degree factors. 

More precisely, algorithm 4 sorts the list of the unknown occurences iiij ac- 
cording to the increasing jdi. The nullities are then computed until there remain 
fewer than a fixed number T of unknowns to be determined by combinatorial 
search. 



Algorithm 4: Nullity-comb-search 
Data: A: an n x n matrix over a finite field, 
Data: P^. the irreducible factors of P^ in , 
Data: T: a static threshold 

Result: m;: the multiplicities of each Pj in P c ^ ar . 
begin 

£ = {(*> i)A = L..k,j = l...ei} 

Sort £ according to the increasing values of jdi 

while (#£ >T)do 

Pop (i,j) from £ 

Compute v>ij = n — rank(P/ (A)) 

for % = 1 . . . k do 

Let ji be the largest index s.t. Vij i is computed 
if ji < ej then 

|_ Compute Vij i+ \ = n — rank(P/ i+1 (A)) 
for k = 1 . . . ji do Compute using cor. 2.6 
Combinatorial-search^, (d u . . .,d k ), (ji, ■ ■ ■ ,jk)) 
for i = 1 . . . k do m ; = j n i,j 
return m = (mi, . . . , m^) 
end 



The combinatorial search has exponential complexity. The threshold T must 
be small. In experiments, we found that T = 5 was the best choice for various 
matrices. In the case of numerous factors with large degree, this approach is of 
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reduced effectiveness. We propose in the next section how to combine it with a 
third algorithm. 



3.2 Nullity method and system resolution 

The index calculus method also enables the design of a hybrid algorithm. If the 
multiplicities of some factors have already been computed by another method, 
we can limit the system to the unknown multiplicities only, thus reducing its 
dimension. 

Suppose that the multiplicities rrii of the factors Pi for i e C are already 
known, then equation (8) reduces to 



mj log^ (A*)) = log(dct(A,J - A)) 

3<tC 

— rrij log(Pj mod q — 1. 



(9) 



A first simple hybrid approach is the following: the method of the nullity 
(section 2.1) is applied to every degree one factor with multiplicity one in the 
minimal polynomial, and the remaining factors are left to the index calculus 
method. 

This approach is always worthy since the computation of the rank of Pj(A), 
for a degree one polynomial Pi is cheaper than the computation of det(\il — A) 

A second hybrid approach also introduces a combinatorial search to this 
algorithm: the nullity method still handles the t degree one factors as previously. 
The remaining factors Pi are sorted by decreasing degree. For a convenient 
choice of s, a list of every possible assignment for the multiplicities of the first 
s factors is determined, using a combinatorial search. Then for each partial 
assignment, the multiplicities of the remaining factors are determined by the 
resolution of an index calculus system of the form: 



mj log(P,-(Ai)) = log g (det(A l J - A)) 

i 

log g (j[Pp (Xi)j modg-lVi. 



(10) 



For each partial assignment, the system resolutions share the same matrix 
B. Therefore the expensive part of it, namely the Gaussian elimination, can be 
performed only once at cost 0(m 3 ), where m = k — s — t. There only remains to 
solve two triangular systems (in 0(m 2 )) for each possible assignment. Lastly, 
the assignments will be discriminated against each other by a test on the total 
degree. To sum up, this techniques makes it possible to balance the cost of the 
computations of determinants, and the cost of the system solving, by reducing 
the dimension of the system, but increasing the dimension of its right hand side. 
The most appropriate value for s has to be determined dynamically, according 
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to the number of possible assignments induced, and using an estimate of the 
cost function of this algorithm: e.g. 

2 

2mnSl + -m 3 + 4m 2 r s 

O 

where t s denote the number of possible assignments for a chosen subset of s 
factors. 

3.3 Index calculus and k"' invariant 

The best known black-box algorithm to compute the Frobcnius normal form over 
a field is given by Villard in [28] . It is proved that computing the fcth invariant 
factor of a matrix reduces to the computation of a minimal polynomial of the 
input matrix with a rank k additive perturbation. Using a binary search tech- 
nique, the algorithm only performs /ilog(n) such computations, where /i is the 
number of distinct invariant factors of the matrix. Since /i is smaller than 3y / n/2 
and an invariant factor can be recovered using 0(n) matrix vector products, this 
corresponds to a total number of 0(n 3 / 2 log(n)) matrix-vector products and an 
additional cost of 0(n 5 / 2 log 2 (n)loglog(n)) arithmetic operations. 

We propose in algorithm 5 an alternative approach combining the index 
calculus method with computations of individual invariant factors. 



Algorithm 5: black-box-charpoly 
Data: A: an n x n matrix over a finite field K 
Result: P^ al or "fail" 
begin 

/i = InvFact(l) 

Factor /i = JJ^ =1 ff * using Cantor-Zassenhaus 
Set S = {P 1 ,...,P k } andj = 2 
while (#5 > VE) do 
fj = InvFact(j) 
forall P, e S do 

Compute a s.t. gcd(ff , fj) = P t a 
if a = then S = S\{Pi} 
|_ mi += a 

return Index- calculus(A, (Pi),S, Yljgs Pj" 3 ) 
end 



The idea is to reduce the dimension of the index calculus system to y/n by 
computing a few of the first invariant factors of the matrix. 

After each computation of an invariant factor $, the multiplicity of each 
irreducible polynomial Pi is updated, and those Pi that are no longer in $ are 
removed from the list of the factors with unknown multiplicity. 
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The while loop is executed at most y/n times. Otherwise, there would be 
more than \fn invariant factors having more than yfn irreducible factors, and 
the total degree would be larger than n. 

Now the condition of exit for this loop ensures that the order of the linear 
system will be smaller than ^fn. Therefore only sjn determinants will be com- 
puted and the overall number of blackbox matrix- vector products is 0(n^/n) 

The remaining multiplicities are then determined by the index calculus 
method described previously, requiring at most 0(n 3 / 2 Q) applications of the 
matrix to a vector. 

Under the conditions of validity for the index calculus algorithm, this heuris- 
tic improves on the computation time of Villard's algorithm by a logarithmic 
factor. 

4 Lifting over the integers 

Storjohann gives in [26] a method for the computation of the Frobenius normal 
form of a black-box integer matrix. It is based on a computation of the minimal 
polynomial over Z and on a computation of the Frobenius normal form over 
a prime field. Then a gcd-free basis for the invariant factors over Z p [X] is 
computed and lifted over 1\X\. 

We use the same idea but just for the characteristic polynomial, and not for 
all the invariant factors. It is thus simpler since we don't need to ensure that the 
Frobenius form of A modulo p equals the integer Frobenius form reduced modulo 
p. We just need ensure that the minimal and characteristic polynomials of A 
modulo p equal the minimal and characteristic polynomials over the integers 
reduced modulo p. 

The goal of the following algorithm 6 is to compute the integer characteristic 
polynomial from the integer minimal polynomial and the characteristic polyno- 
mial modulo some prime p, obtained via the previous sections. Algorithm 6 is 
just a simplification of that of [26]: 

The integer minimal polynomial is computed via [9, Theorem 3.3] with an 
O(sdCt) probabilistic complexity, where s is the size of its integer coefficients and 
d, its degree. The characteristic polynomial modulo p is computed via algorithm 

5 with an 0(n 15 f2) complexity and (D~(n 2 - 5 ) extra field operations. Then the 
squarefree part [17] and the Henscl lifting of the gcd-free basis takes (D~(nk) word 
operations with fast integer and polynomial arithmetic [16, Theorem 15.18]. 

The size of the coefficients of the integer minimal polynomial is bounded 
in the worst case by s < §(log 2 (n) + log 2 (| \A\ | 2 ) + 0.212) where \\A\\ is the 
largest entry in absolute value of the matrix A, see [7, Lemma 2.1]. When 
the matrix entries are of constant size, s — (D~(n) and as the degree of the 
minimal polynomial is bounded by n, the dominant asymptotic cost is that of 
the integer minimal polynomial computation. This result is already in [26]. In 
practice, however, the coefficients of the minimal polynomial are often much 
smaller than the bound and than those of the integer characteristic polynomial. 
Furthermore, the degree of the minimal polynomial can be extremely small, 
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Algorithm 6: Gcd-free lifting of the characteristic 

polynomial 

Data: A: an n x n integer matrix, P^ in its minimal polynomial over Z, 
Data: p: a prime number, 

Data: P^ the characteristic polynomial of A mod p, 
Result: P^ ml the characteristic polynomial of A over Z 
begin 

Compute S the squarefree part of P^in 

S = S mod p, PA n = PA n mod p 

Compute a modular gcd-free basis (gi, ... , g x ) of (S, P^ in , -P c h ar ), 
together with exponents (/ii, . . . , /i x ) such that P c ^ ar = J| g^ 
Apply Hcnscl lifting on the basis (g~i, . . . ,g x ) to produce {g\, . . . ,g x ) 
so that S = g 1 ...g x mod p k and (g u . . . , g x ) = (g u . . . , g x ) mod p 
return P^ MT = Y[g^ 
end 



especially for structured or sparse matrices (see e.g. homology matrices in [9]). 
In those cases, the dominant cost will be the computation of the characteristic 
polynomial modulo p. Then, our algorithm enables faster computations since a 
factor of ^/n has been gained, as illustrated in table 1. 

A supplemental constraint can be introduced by computing dct( XI — A), at 
random integer A. Using e.g. [10, Theorem 4.2] with [13], 0{s/n) of these can be 
done to speed-up the modular adaptive search, without increased complexity. 

5 Experimental comparisons and applications 

We have implemented some of the algorithms presented in the previous sections 
using the LinBox 2 library for the black-box computation of minimal polynomi- 
als, ranks and determinants. In a first approach, we replaced the computation 
of the gcd-free basis of algorithm 6 by a factorization into irreducible factors, 
using Hensel lifting. This algorithm is more expensive in the worst case, but 
the efficient implementation by NTL 3 makes it practicable in numerous cases. 

This work was partly motivated by an application from graph theory. For 
a graph X on n vertices with vertex set V(X) and edge set E(X) the k-th 
symmetric power X k is the graph with the (^) fc-subsets of V(X) as vertices 
and with two such fc-subsets adjacent if their symmetric difference is in E(X). 

Graph theorists are interested in the spectrum of such graphs (defined as 
the spectrum of their adjacency matrix) since they are closely related to the 
description of their isomorphism class. More precisely, if a certain power fc is 
found, such that the symmetric fcth power of a graph describes its isomorphism 

2 www. linalg.org 
3 www . shoup . net /nt 1 
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class, this would provide a polynomial algorithm to solve the graph isomorphism 
problem. 

This motivated the team of Audenaert, Godsil, Royle and Rudolph to study 
in [1] the spectrum of symmetric powers of a class of graphs: the strongly regular 
graphs. They prove that there exist infinitely many graphs having co-spectral 
symmetric squares. But concerning the symmetric cubes, no pair of graph is 
known to have co-spectral symmetric cubes until now. 

We helped Royle to investigate further the computation of the characteristic 
polynomial of the symmetric cubes of strongly regular graphs. He was able to 
test the first 72 cases corresponding to the graphs with fewer than 29 vertices. 
Using our implementations available in LinBox, we have been able to test the 
36 582 graphs with fewer than 36 vertices and check that there is no pair of 
graphs among them having cospectral symmetric cubes. 

We used the matrices of this application to benchmark the implementations 
of the previously presented algorithms. These matrices are sparse and sym- 
metric, and therefore especially suited to black-box computations. Moreover, 
several parameters such as the degrees of their minimal polynomials or the av- 
erage number u of nonzero elements per row vary among the matrices. The 
matrices EX1 , EX3 , EX5 correspond respectively to the symmetric cubes of the 
strongly regular graphs with parameters (16,6,2,2), (26,10,3,4) and (35,16,6,8). 
Their dimensions are respectively 560 = (g 6 ) , 2600 = ( 2 3 6 ) and 6545 = ( 3 3 5 ) . The 
matrices EX2 and EX4 correspond to different graphs but with similar parameters 
as EX1 and EX3. 

All the matrices used in the experiments, including adjacency matrices of the 
symmetric powers, are avaible on-line in the Sparse Integer Matrices Collection 1 . 
In particular we used, in the following tables, matrices from the SPG, Forest, 
Trefethen and Homology sections of the collection. When the tested matrix 
was not square, we considered the square matrix obtained by padding it with 
zeroes. 

We report in table 1 the computation time of the different modules described 
in this paper. For each matrix, two computations are compared: they share the 
computation of the minimal polynomial over Z. Then the determination of the 
multiplicities is done either by the combination of the nullity algorithm and the 
combinatorial search (with the threshold T set to 5), or by the index calculus 
method. 

We first note that the determination of the multiplicities may be the domi- 
nant operation when the degree of the minimal polynomial is small, as for the 
matrix EX1. This makes the motivation for this study obvious. For this task 
either method, nullity or resolution of the logarithmic system, can be the most 
competitive option, depending on the structure of the irreducible factors. This 
advocates for the adaptive approach of algorithm 5 combining both methods, 
and the computation of the fcth invariant factor. 

In order to emphasize the improvement of the black-box determination of the 
multiplicities over dense methods, we now compare it to the alternative tech- 

4 ljk. imag.fr/membres/Jean-Guillaume.Dumas/SIMC 
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Matrix 


EX1 


EX2 


EX3 


EX4 


EX5 


n: dimension 


560 


560 


2600 


2600 


6545 


d: deg (P min ) 


54 


103 


1036 


1552 


2874 


u>: sparsity 


15.6 


15.6 


27.6 


27.6 


45.2 


Z-Minpoly 


0.11s 


0.26s 


117s 


260s 


5002s 


Z[X] factorize 


0.02s 


0.07s 


9.4 


18.15 


74.09s 


Nullity /comb. 


3.37s 


5.33s 


33.2s 


30.15s 


289s 


Total 


3.51s 


5.66s 


159.4s 


308.1s 


5366s 


Index calc. 


3.46s 


4.31s 


64.0s 


57.0s 


647s 


Total 


3.59s 


4.64s 


190.4s 


336.4s 


5641s 



Table 1: Computation time for tasks of the integer adaptive algorithm on a 
Pcntmm4 (3.2 GHz; 1 Gb) 

nique presented in [8, §4.2.2]. This also relies on the computation of the minimal 
polynomial in 7L\X\ and its decomposition into irreducible factors. But the mul- 
tiplicities are then obtained using one dense computation of the characteristic 
polynomial in a randomly chosen finite field. It is therefore not anymore a black- 
box algorithm. We will denote it by dchar. Comb is the nullity-combinatorial 
search algorithm, and ind is the index calculus method. A denotes 08blocks, 
B is ch5-5.b3, and T is Tref500 from the Sparse Integer Matrices Collection. 



Matrix 


n 




dchar 


null- comb 


ind . 


A 


300 


1.9 


0.32s 


0.08s 


0.07s 


AA T 


300 


2.95 


0.81s 


0.12s 


0.12s 


B 


600 


4 


4.4s 


1.52s 


1.97s 


BB T 


600 


13 


2.15s 


3.96 


7.48s 


TF12 


552 


7.6 


6.8s 


5.53s 


5.75s 


mk9b3 


1260 


3 


31.25s 


10.51s 


177s 


Tref500 


500 


16.9 


65.14s 


25.14s 


25.17s 



Table 2: Integer black-box approach for multiplicities on an Athlon (1.8 GHz; 
2 Gb) 

Table 2 shows the improvement of the black-box approach for several ma- 
trices coming from different applications. Once again, the structure of the ir- 
reducible factors of the minimal polynomials cause various behaviors for each 
variant. For example the times of Index-calculus are similar to those of 
Nullity-comb-search, sometimes better but also sometimes much slower, as 
for the matrices BB T and mk9b3. 
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6 Conclusion 



We developed several ways to recover the multiplicities of the factors of the char- 
acteristic polynomial from a factorization of the minimal polynomial. Over a 
finite field hybrid heuristics arc proposed, that compete with the best theoretical 
complexity. Over the ring of integers, our approach enables fast computations 
particularly when the coefficients or degree of the minimal polynomial arc small. 
This is illustrated on a family of strongly regular graphs, in order to verify that 
there are no symmetric co-spectral cubes. 

Further studies on the theoretical complexity remain to be done, and could 
lead to better implementations in practice. In particular, a recent algorithm for 
dense matrices [24] might be adapted for black-box matrices. In this regard, 
extending the block projections of [14] to the case of similarity transformations 
would play a crucial role. 
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